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Abstract 

Prey predator population dynamics has been frequently 
studied in literature to model many robust but contrastive 
biological populations. In this paper, based on prey predator 
population dynamics, a nonlinear model to predict the 
behavior of virus treatment of cancer cells is suggested. Based 
on the goals of virus curing and consideration of some 
practical issues, the suggested model has the capability to 
cover these following cases: both cancer cells and viruses 
vanish, one is eliminated by the other, or both survive in a 
robust periodically convergent manner or both become 
divergent. 

Relative to the values of some parameters such as the virus 
intrinsic growth rate, carrying capacity, capturing rate, half 
saturation constant, maximal tumor cells growth rate, and 
tumor cells mortality rate, and the initial population values of 
both tumor cells and viruses, some examples are solved. These 
examples exhibit logical and feasible behavior for the system. 
They show under which conditions the virus therapy fails or 
succeeds; as well when the bifurcation occurs or when the 
trajectories converge to a limit cycle. 
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Introduction 

Cancer is one of biggest causes of death, thus finding 
ways to beat cancer is really vital. Many researches can 
be found in literature concerning about the cancer 
treatment, tumor cell behaviors, clinical cares, etc; but 
because the main source of change from cells to tumor 
and defected ones is not clearly known, the treatment of 
cancer is almost based on the reduction of the 
destructive effects of cell behaviors. Virus therapy is an 



almost new method used for some types of cancer to 
destroy tumor cells by injecting killer viruses into the 
core of infected cells. Practically, this method is based on 
a trial and error process whereas a model is constructed 
to investigate the behaviors of tumor cells, which is 
really important to improve the researches in this area 
and to help the medicos to treat such diseases. By means 
of the model, the cell behaviors can be investigated and 
then the treatment process gets controlled. 

In the past decades, some researchers have tried to 
utilize the virus therapy method to cancer treatment. 
However, due to the lack of molecular biology 
statements for the viruses, until 1990 these techniques 
didn't attract sufficient attention. Based on recent 
accomplishments in molecular biology of the viruses, the 
conventional methods in cancer treatment such as 
chemotherapy and radiotherapy can be replaced by 
virus therapy. The main superiority of the last method 
with respect to chemotherapy and radiotherapy is that 
the former kills the defected cells only see for example 
(Daniel Cervantes et. al. 2008) 

Many mathematical models have been suggested for 
virus therapies of cancer tumor; see for example (Artem. 
S. et al. 2006, Dominik Wodarz, 2001. Zeljko Bajzer, et.al, 
2008. Faina S, et al, 2007. Georgy P Karev, et al, 2006). 
With respect to their advantages, they cannot cover the 
whole dynamics of the process because in some of them, 
the death of all cancer cells cannot be reached or, for 
example, the coexistence of cancer cells and viruses is 
not probable. Although the reference (Dominik Wodarz. 
et al. 2009), presented an almost good model but all the 
probable coexistence cases of tumor cells and viruses 
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were not covered. In reference (Manju Agarwal et. al. 
2011), a good coverage over the probable cases is also 
seen. SH Thorne and CH Contag (2008) studied the 
improvement of therapeutic benefits of cell-based 
delivery when the biological characteristics of oncolytic 
viruses and immune cells are integrated. Another 
approach studying the coexistence of two predators on 
one prey (see Irakli Loladze et. al. 2004) may be used to 
model the virus therapy of cancers with two types of 
viruses. 

The Modified Prey Predator Model 

Based on the famous prey predator model that was 
presented in (Bingtuan LI and Yang Kuang, 2007), a 
modified model is proposed in this article to predict the 
dynamics of the population of both viruses and tumor 
cells. Here the following equations are in use to describe 
the dynamics: 
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Where x and y are the population of prey (tumor cells) 
and predator (virus), respectively. In addition, r, k, c, m, f, 
and d , the positive constants demonstrate virus intrinsic 
growth rate, carrying capacity, capturing rate, half 
saturation constant, maximal tumor cells growth rate, 
and tumor cells mortality rate, respectively. For 
simplicity, using following non dimensional parameters, 

x —> Kx, y —> Ky I mj —> mt I c (3) 
one can reach such the following format: 
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If both x and y take positive values, relations (4) and (5) 
are reduced to: 
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x\t) = ax(l - x)(x + y) - xy 

y'(t) = -Py(x+y) + Kxy 
Using a change of variable, 

t -^>t(x+ y) , 
and then Briot-Bouquet's transformation, 

x^x, y—>yx,t—>tlx, (10) 
one can reach: 

x\t) = x[a - ax - (1 - a)y - axy] (11) 
y\t) = y[(K-a- P) + ax + 
(I- a- j3) + axy] 
Replacing x -> x/a in (11) and (12), one can obtain: 

x'(t) = x[a - x - (1 - a)y - xy] (13) 
y\t) = y[(K-a-P) + x + 

(l-a-j3) + xy] 
a and v = k - a - [3 (or equivalently a and k) are utilized 
as our unfolding parameters while [3 is fixed. In (13 & 14) 
there are two second order terms whose coefficients 
depend on a. Then these terms are decomposed and 
rewritten (13 & 14) as: 

x{t) - x(a — x — y) + x(ay - xy) (15) 

y'(t) = y(y + x + (l-J3)y) (16) 
Using the transformations 

x gx , y -^sy , a = sv x , 
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and rescaling time t -* tff, we convert (15 & 16) to: 

x'(t) = x[v 1 -x-y] + £(v 2 xy-x 2 y) (18) 



y'd) = y[ ^ V, + x + (1 - /?)] + 

s(v 2 y-v l y 2 + y 2 x) 



(19) 
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The system of equations (18 & 19) describes the 
nonlinear model of the viruses and cancer cell 
interaction in this work. 

Examples 

The parameters v L , v 2 , and 8 in (18 & 19) related to r, k, 
c, m, f, and d have been introduced previously as virus 
intrinsic growth rate, carrying capacity, capturing rate, 
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half saturation constant, maximal tumor cells growth 
rate, and tumor cells mortality rate, respectively. Here, 
based on these parameters, the stationary points and the 
behavior of the system are obtained and drawn in phase 
plane. 

Case I The Virus Therapy Fails 

By means of the values presented in Table 1, and the 
solution to the equations, the behavior drawn in Fig. 1 is 
confronted. 

TABLE 1 THE VALUES OF PARAMETERS USED IN CASE 1 
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FIG. 1 THE BEHAVIOR OF THE SYSTEM IN PHASE PLANE FOR 
CASE I 

( x: population of tumor cells; y: population of viruses) 

According to this figure, only two saddle points 
constructed at (0, 0) and (40, 0) are the stationary points 
of the system. Fig. 1 shows also that the virus therapy in 
this case fails such that the trajectories initiated from 
each initial condition (all points in the phase plane 
except y axis) terminate to stationary point (40,0), 
meanings that after sufficient amount of time, the 
viruses are vanished by the tumor cells. 

Case II The Virus Therapy Succeeds 

By means of the values presented in Table 2, as well as 



the solution to the equations, the behavior drawn in Fig. 
2 is confronted as well which shows a bifurcation. In 
other words, since the factor p (see (6) also) reduces to 
smaller values, the number of the stationary points and 
the whole behavior completely change. 

TABLE 2 THE VALUES OF PARAMETERS USED IN CASE 2 
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FIG. 2 THE BEHAVIOR OF THE SYSTEM IN PHASE PLANE FOR 
CASE II 

( x: population of tumor cells; y: population of viruses) 

Based on this figure, a source at (0, 0), a saddle point at 
(40, 0) and a sink at (0,159.8721) can be seen, as well as 
shows that the virus therapy in this case succeeds such 
that the trajectories initiated from each initial condition 
(all points in the phase plane except ones on the 
convergent manifold of the saddle point) terminate to 
stationary point (0,159.8721), meaning that after 
sufficient amount of time, the tumor cells are vanished 
by viruses. 

Case III Coexistence of Both Viruses and Tumor Cells 

By means of the values presented in Table 3, and the 
solution to the equations, the behavior drawn in Fig. 3 is 
confronted as well which shows another bifurcation also 
the famous robust limit cycle usually seen in prey 
predator models. In other words, since the factor p (see 
(6) also) reduces to smaller values but v 2 becomes larger, 
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the number of the stationary points and the whole 
behavior of system completely change. 

TABLE 3 THE VALUES OF PARAMETERS USED IN CASE 3 
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FIG. 3 THE BEHAVIOR OF THE SYSTEM IN PHASE PLANE FOR 
CASE III 

(x: population of tumor cells; y: population of viruses) 

Based on this figure, three saddle points at (0, 0), (40,0) 
and (0,60) and a stable focus at (12.062, 28.4) can be seen. 

Fig. 3 shows also that the virus therapy in this case 
results in the coexistence between the viruses and the 
tumor cells. In this case the trajectories initiated from a 
point lain in the triangle constructed by three saddle 
points converge to the stable focus that presents a 
chronic illness. For other initial conditions, the tumor 
cells are reduced slowly but not completely removed. 

Quantitatively, changing the parameters does not lead to 
significant changes in behavior, whereas the rate of 
response, stationary points, areas of attraction and 
region of the robustness of limit cycle will change. In Fig. 
4, a slower dynamics and in Fig.5 the evolution of the 
limit cycle for the case 3 are drawn 



15 20 25 30 35 40 45 



( x: population of tumor cells; y: population of viruses) 

It is noted that the area constructed in the triangle in Fig. 
4 is in direct relation to the squared value of v z . This may 
be utilized in validation of the model with experimental 
studies and prediction on the behavior for some 
standard practical cases. 
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FIG. 5 THE EVOLUTION OF THE LIMIT CYCLE OF FIG. 4 
( x: population of tumor cells; y: population of viruses) 

Local Stability Analysis of Equilibrium Points 

Local stability analysis is concerned with the stability of 
a nonlinear system near its equilibrium points. It is a 
formalization of the intuition that a nonlinear system 
should behave similarly to its linearized approximation 
for small range motions. See (Jean-Jacques E. Slotine and 
Weiping LI. 2001, Erwin Kreyszig, 2010) for more details. 

Here for our model, to check the stability near 
equilibrium points, the variational matrix is constructed 
by differentiating the right hand side of Equations 
(18&19) with respect to , x and y. The entries of the 
general variational matrix V for point P become: 



V(p) = 



A B 
C D 



Where 



FIG. 4 THE BEHAVIOR OF THE SYSTEM IN PHASE PLANE 
SIMILAR TO CASE III BUT WITH SMALLER p AND LARGER 



A = v l -2x-(sv l -2£K-l)y 

B = sv x x — sx 2 — X 
C = y + ey 
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(21) 
(22) 
(23) 
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D= 2(1 P\ .+x+2(l-P)y + ev 2 

2-fi (24) 

-2sv x y + 2sxy 
Now for studied cases, the local stabilities in equilibrium 
points are investigated. To explore the local stability of 
trivial equilibrium points, the eigenvalues of V (P) at 
stationary points are computed and analyzed. 



Case I 

1) Local Stability Analysis of P(0,0) 
The V (P) and its eigenvalues are: 

"40 



V(p(40,0)) = 



-240 



X = 40 , X = -240 
So we have a saddle point at P (0, 0) 

2) Local Stability Analysis of P(40,0) 

The V (P) and its eigenvalues are: 



V(p(40,0)) = 



-40 -40 
-200 
X = -40 , X = -200 

So we have a nodal sink point at P (40, 0) 
Case II 

1) Local Stability Analysis of P(0,0) 
The V (P) and its eigenvalues are: 



V(p(0,0)) = 



40 
80 



A = 40 , A = 80 
So we have a nodal source point at P (0, 0) 

2) Local Stability Analysis of P(0,159.8721) 

The V (P) and its eigenvalues are: 

"- 119.8082 - 9.9057 e- 1 

160.1277 - 80 

A = -119.8082 , A = -80 
So we have a nodal sink point at P (0,159.8721) 

3) Local Stability Analysis of P(40,0) 
The V (P) and its eigenvalues are: 



V(p(0, 159.8721)) = 



-40 -40 
■5.0741^-7 120 

A = -40 , A = 120 
So we have a saddle point at P (40,0) 

Case III 

1) Local Stability Analysis of P(0,0) 
The Eigenvalues of V (P) at point P(0,0) is given by: 

"40 

V(p(0,0)) = 

[0 -22.8589_ 

A = 40 , A = -22.8589 
So we have a saddle point at P (0, 0) 

2) Local Stability Analysis of P(11.8307,28.5718) 
The V (P) and its eigenvalues are: 

V(p(l 1.8307,28 .5718)) = 
"-11.9997 -11.664" 
28.9799 11.0263 _ 
2 = -0.486697 +14.3344/ 
So we have a spiral sink or stable focus at this point. 

3) Local Stability Analysis of P(0,60.1498) 
The V (P) and its eigenvalues are: 

I"- 18.9462 

V(p(0,60.1498)) = 

[ 61.9589 22.8569 
X = -18.9462 , A = 22.8569 
So we have a saddle point at P (0, 60.1498). 

4) Local Stability Analysis of P(40,0) 
The V (P) and its eigenvalues are: 

-40 40 
•3.5752^-7 17.1431 

A = -40 , A = 17.1431 
So we have a saddle point at P (40, 0). 

Conclusion 

In this paper, a modified model based on the prey 
predator population dynamics is presented with the 
capability to predict and investigate the practical clinical 
cases. Relative values of some parameters to the virus 



v(pmoy) = 
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intrinsic growth rate, carrying capacity, capturing rate, 
half saturation constant, maximal tumor cells growth 
rate, and tumor cells mortality rate, and the initial 
population values of both tumor cells and viruses, some 
examples are solved which exhibit logical and feasible 
behavior for the system. They show under which 
conditions the virus therapy fails or succeeds; also when 
the bifurcation occurs or when the trajectories converge 
to a limit cycle. 
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